#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* THE GLORIOUS TRINITY GEOMETRY LIBRARY version 1.1 */

/* In general, a return code 2 means strictly within; 1 on the boundary;
 * 0 outside.
 * Exactly one preprocessor statement USE_FLOAT or USE_RATIO should be defined.
 * Alternatively, edit these yourself.
 * Methods requiring floating point are 'angle_bisector', 'incentre',
 * 'intersect_*' & 'tangent_*'.
 * At infinity, anything might happen, unless you only use 'cross' and
 * 'full_dot'.
 * Reverse dependencies assume you have typed up ALL the Core
*/

// Possible options
// #define USE_FLOAT
// #define USE_RATIO

/* THE CORE */

// md5begin core
#ifdef USE_RATIO
  typedef long long ptype;
#endif
#ifdef USE_FLOAT
  typedef double ptype;
#define EPSILON 1.0e-15
#endif

typedef struct { ptype x, y, z; } Point;
typedef Point Line;
// A rational is a Point with y = 0, representing x / z
typedef Point Rational;

ptype ptypeabs(ptype x) {
  return x >= 0 ? x : -x;
}

#ifdef USE_RATIO
ptype ptypegcd(ptype x, ptype y)
{
  ptype t;
  x = ptypeabs(x);
  y = ptypeabs(y);
  while(y > 0)
  {
    t = x % y;
    x = y;
    y = t;
  }
  return x;
}
#endif

Point mk_point(ptype x, ptype y, ptype z) {
  Point t;
  t.x = x, t.y = y, t.z = z;
  return t;
}

Point reduce(const Point &t) {
  if(t.x == 0 && t.y == 0 && t.z == 0)
    return t;
#ifdef USE_RATIO
  ptype r = ptypegcd(ptypegcd(t.x, t.y), t.z);
#endif
#ifdef USE_FLOAT
  ptype r = (t.z == 0 || ptypeabs(t.x / t.z) > 1.0 / EPSILON) ?
      ptypeabs(t.x) + ptypeabs(t.y) : ptypeabs(t.z);
#endif
  if(t.z < 0)
    r = -r;
  return mk_point(t.x/r, t.y/r, t.z/r);
}
// md5end

/* ARITHMETIC */

// md5begin cross
Point cross(const Point &t, const Point &u) {
  return reduce(mk_point(t.y * u.z - t.z * u.y,
                         t.z * u.x - t.x * u.z,
                         t.x * u.y - t.y * u.x));
}
// md5end

// md5begin full_dot
ptype full_dot(const Point &t, const Point &u) {
  return t.x * u.x + t.y * u.y + t.z * u.z;
}
// md5end

// md5begin vec_add
Point vec_add(const Point &t, const Point &u) {
  return reduce(mk_point(t.x * u.z + u.x * t.z,
                         t.y * u.z + u.y * t.z,
                         t.z * u.z));
}
// md5end

// md5begin vec_sub
Point vec_sub(const Point &t, const Point &u) {
  return reduce(mk_point(t.x * u.z - u.x * t.z,
        t.y * u.z - u.y * t.z,
                         t.z * u.z));
}
// md5end

// md5begin scalar_mult
Point scalar_mult(const Point &t, const Rational &k) {
  return reduce(mk_point(t.x * k.x, t.y * k.x, t.z * k.z));
}
// md5end

// md5begin cart_dot
Rational cart_dot(const Point &t, const Point &u) {
  return reduce(mk_point(t.x * u.x + t.y * u.y, 0, t.z * u.z));
}
// md5end

// md5begin norm
Rational norm(const Point &t) {
// Requires cart_dot
  return cart_dot(t, t);
}
// md5end

// md5begin cart_cross
Rational cart_cross(const Point &t, const Point &u) {
  return reduce(mk_point(t.x * u.y - t.y * u.x, 0, t.z * u.z));
}
// md5end

// md5begin are_equal
char are_equal(const Point &t, const Point &u) {
// Note values should have been reduced recently to call this
#ifdef USE_RATIO
  return (t.x == u.x && t.y == u.y && t.z == u.z)
      || (t.x == -u.x && t.y == -u.y && t.z == -u.z);
#endif
#ifdef USE_FLOAT
  return (ptypeabs(t.x - u.x) < EPSILON)
      && (ptypeabs(t.y - u.y) < EPSILON)
      && (ptypeabs(t.z - u.z) < EPSILON);
#endif
}
// md5end

// md5begin compare_ratio
ptype compare_ratio(const Rational &t, const Rational &u) {
  return t.z == 0 ? u.z == 0 ? 0 : 1 : u.z == 0 ? -1 : (t.x * u.z - u.x * t.z);
}
// md5end

// md5begin side_of
char side_of(const Point &x1, const Point &x2, const Point &y) {
// Requires vec_sub, cart_cross
// LEFT = -1, RIGHT = 1, ON = 0
  Point x2b = vec_sub(x2, x1);
  Point yb = vec_sub(y, x1);
  Rational p = cart_cross(yb, x2b);
  return p.x > 0 ? 1 : p.x < 0 ? -1 : 0;
}
// md5end

/* CONTAINMENT */

// md5begin containment
typedef enum { DontX = 0, X2Ends, X1End, XNormal, XInf } XType;
char seg_contains(const Point &x1, const Point &x2, const Point &y)
// Requires vec_sub, cart_dot
{
  Point s = cart_dot(vec_sub(x1, y), vec_sub(x2, y));
  if(s.x == 0)
    return 1;
  if(s.x*s.z > 0)
    return 0;
  return 2;
}

XType segmentX(const Point &x1, const Point &x2, const Point &y1,
               const Point &y2, Point *z)
// Requires (local) cross, are_equal, seg_contains; (inherited) vec_sub,
// cart_dot
// Not as bad as it looks, especially if the VERY special case of collinearity
// cannot occur; removes are_equal dep.
{
  Line lx = cross(x1, x2), ly = cross(y1, y2);
  Point i;
  char s, s1, s2;

  if(are_equal(lx, ly))
  // Skip this if segments not on same line --
  {
    s = seg_contains(y1, y2, x1) | seg_contains(y1, y2, x2)
        | seg_contains(x1, x2, y1) | seg_contains(x1, x2, y2);
    if(s & 2)
      return XInf;
    if((are_equal(x1, y1) && are_equal(x2, y2))
         || (are_equal(x1, y2) && are_equal(x2, y1)))
       return XInf;
    if(s & 1)
      return X2Ends;
    return DontX;
  }
  // -- to here

  i = cross(lx, ly);
  s1 = seg_contains(x1, x2, i), s2 = seg_contains(y1, y2, i);

  if(s1 == 0 || s2 == 0)
    return DontX;
  if(z != NULL)
    *z = i;
  if(s1+s2 == 3)
    return X1End;
  if(s1 == 1 && s2 == 1)
    return X2Ends;
  return XNormal;
}
// md5end

// md5begin poly_contains
char poly_contains(Point r, Point *xs, int n)
// Requires (local) seg_contains, side_of, full_dot; (inherited) vec_sub,
// cart_dot, cart_cross
// Uses winding #s
{
  int ii, wn = 0;
  char s;
  Line l = mk_point(0, r.z, -r.y);
  for(ii = 0; ii < n; ii++)
  {
    s = side_of(r, xs[ii], xs[(ii + 1) % n]);
    if (s == 0 && seg_contains(xs[ii], xs[(ii + 1) % n], r))
      return 1;
    // On the poly
    if ((full_dot(xs[ii], l) >= 0) == (full_dot(xs[(ii + 1) % n], l) < 0))
      wn += s;
  }
  return wn == 0 ? 0 : 2;
}
// md5end

/* METRICS */

// md5begin tri_area
Rational tri_area(const Point &x, const Point &y, const Point &z)
// Requires cart_cross, vec_sub
// NOTE signed area
{
  Point a = cart_cross(vec_sub(y, x), vec_sub(z, x));
  a.z *= 2;
  return reduce(a);
}
// md5end

// md5begin dist2_pt_line
Rational dist2_pt_line(const Point &x, const Point &l)
// Requires full_dot
{
  Point t;
  t.x = full_dot(x, l);
  t.y = 0;
  t.x *= t.x;
  t.z = x.z*x.z*(l.x * l.x + l.y * l.y);

  return reduce(t);
}
// md5end

// md5begin dist2_pt_sgmt
Rational dist2_pt_sgmt(const Point &x1, const Point &x2, const Point &y)
// Requires (local) cart_dot, vec_sub, norm, dist2_pt_line; (inherited)
// full_dot, cart_dot
{
  Point l = cross(x1, x2);
  Point d1 = cart_dot(vec_sub(x2, x1), vec_sub(y, x1));
  Point d2 = cart_dot(vec_sub(x1, x2), vec_sub(y, x2));

  if(d1.x < 0)
    return norm(vec_sub(y, x1));
  if(d2.x < 0)
    return norm(vec_sub(y, x2));

  return dist2_pt_line(y, l);
}
// md5end

// md5begin dst2_sgmt_line
Rational dist2_sgmt_line(const Point &x1, const Point &x2, const Line &l)
// Requires full_dot, dist2_pt_line, compare_ratio
{
  if(full_dot(x1, l) * full_dot(x2, l) <= 0)
    return mk_point(0, 0, 1);
  Point d1 = dist2_pt_line(x1, l), d2 = dist2_pt_line(x2, l);

  return (compare_ratio(d1, d2) > 0 ? d1 : d2);
}
// md5end

// md5begin dist2_sgmt_sgmt
Rational dist2_sgmt_sgmt(const Point &x1, const Point &x2, const Point &y1,
                         const Point &y2)
// Requires (local) segmentX, dist2_pt_sgmt, compare_ratio; (inherited)
// dist2_pt_line, full_dot, cross, seg_contains, are_equal, norm, vec_sub,
// cart_dot
{
  if(segmentX(x1, x2, y1, y2, NULL) != DontX)
    return mk_point(0, 0, 1);

  Point min = mk_point(1, 0, 0), t;

  t = dist2_pt_sgmt(x1, x2, y1);
  if(compare_ratio(t, min) < 0)
    min = t;
  t = dist2_pt_sgmt(x1, x2, y2);
  if(compare_ratio(t, min) < 0)
    min = t;
  t = dist2_pt_sgmt(y1, y2, x1);
  if(compare_ratio(t, min) < 0)
    min = t;
  t = dist2_pt_sgmt(y1, y2, x2);
  if(compare_ratio(t, min) < 0)
    min = t;

  return min;
}
// md5end

// md5begin poly_area
Rational poly_area(Point *xs, int no_points)
// Requires (local) vec_add, tri_area; (inherited) cart_cross, vec_sub
// NOTE signed area
{
  int ii;
  Point r = mk_point(0, 0, 1), s = r;
  for(ii = 0; ii < no_points; ii++)
    s = vec_add(s, tri_area(r, xs[ii], xs[(ii + 1) % no_points]));
  return s;
}

// md5end

/* CONSTRUCTION */

// md5begin angle_bisector
Line angle_bisector(Point a, Point b, Point c)
// Requires (local) vec_sub, norm, cross, vec_add; (inherited) cart_dot
// NOTE bisector is the one between BA and BC internally
{
  ptype ra, rc;
  a = vec_sub(a, b);
  c = vec_sub(c, b);

  ra = sqrt(norm(a).x);
  rc = sqrt(norm(c).x);

  a.x /= ra;
  a.y /= ra;
  c.x /= rc;
  c.y /= rc;

  return cross(b, vec_add(b, vec_add(a, c)));
}
// md5end

// md5begin incentre
Point incentre(Point a, Point b, Point c)
// Requires (local) angle_bisector, cross; (inherited) vec_sub, norm, vec_add,
// cart_dot
{
  Point l1 = angle_bisector(a, b, c), l2 = angle_bisector(b, c, a);
  return cross(l1, l2);
}
// md5end

// md5begin parallel_line
Point parallel_line(Point a, Point l)
// Requires full_dot
// Through A, || L
{
  l.z = l.z * a.z - full_dot(a, l);
  l.x *= a.z;
  l.y *= a.z;
  return reduce(l);
}
// md5end

// md5begin perp_line
Point perp_line(Point a, Point l)
// Requires (local) parallel_line; (inherited) full_dot
// Line through A, _|_ L
{
  ptype t = l.x;
  l.x = l.y, l.y = -t;
  return parallel_line(a, l);
}
// md5end

// md5begin reflect
Line reflect(Line a, Line b)
// Line b reflected in line a
{
  ptype t = 2 * (a.x * b.x + a.y * b.y);
  ptype u = a.x * a.x + a.y * a.y;

  return reduce(mk_point(b.x * u - a.x * t,
                         b.y * u - a.y * t,
                         b.z * u - a.z * t));
}
// md5end

// md5begin perp_bisector
Point perp_bisector(Point a, Point b)
// Requires (local) perp_line, vec_add; (inherited) parallel_line, full_dot
{
  Point m = vec_add(a, b);
  m.z *= 2;
  return perp_line(m, cross(a, b));
}
// md5end

// md5begin circumencentre
Point circumcentre(Point a, Point b, Point c)
// Requires (local) cross, perp_bisector; (inherited) perp_line, vec_add,
// parallel_line, full_dot
{
  return cross(perp_bisector(a, b), perp_bisector(b, c));
}
// md5end

// md5begin intersect_given_x
char intersect_given_x(Point o, double r, Point x, Point *a1, Point *a2)
// Requires vec_sub, vec_add
{
  double d = (x.x * x.x + x.y * x.y) / (x.z * x.z);
  double t;
  Point y = x;
  if(d > r * r)
    return 0;

  t = y.x;
  y.x = -y.y;
  y.y = t;
  y.z /= sqrt(r * r / d - 1);

  *a1 = vec_add(o, vec_add(x, y));
  *a2 = vec_add(o, vec_sub(x, y));
  return 1;
}
// md5end

// md5begin intersect_circle_circle
char intersect_circle_circle(Point o1, double r1, Point o2, double r2,
                             Point *a1, Point *a2)
// Requires (local) vec_sub, intersect_given_x; (inherited) vec_add
{
  Point p = vec_sub(o2, o1);
  double d = (p.x * p.x + p.y * p.y) / (p.z * p.z);
  if(d >= (r1 + r2) * (r1 + r2) || d <= (r1 - r2) * (r1 - r2))
    return 0;
  p.z *= 2 * d / (d + (r1 * r1 - r2 * r2));

  return intersect_given_x(o1, r1, p, a1, a2);
}
// md5end

// md5begin tangent_circle_point
char tangent_cirle_point(Point o, double r, Point p, Point *a1, Point *a2)
// Requires (local) vec_sub, intersect_given_x; (inherited) vec_add
{
  Point x = vec_sub(p, o);
  double d = ((x.x * x.x + x.y * x.y) / (x.z * x.z));
  x.z *= d;
  return intersect_given_x(o, r, x, a1, a2);
}
// md5end

// md5begin tangent_circle_circle
char tangent_circle_circle(Point o1, double r1, Point o2, double r2, Line *ls)
// Requires (local) vec_add, tangent_cirle_point; (inherited) vec_sub,
// intersect_given_x
// ls is an array large enough for 4 Lines to be put there
{
  Point c, t1, t2;
  t1 = o1;
  t1.z /= r2;
  t2 = o2;
  t2.z /= r1;
  c = vec_add(t1, t2);
  c.z /= (r1 + r2);

  if(!tangent_cirle_point(o2, r2, c, &ls[0], &ls[1]))
    return 0;
  if(r1 >= 0 && r2 >= 0)
    return tangent_circle_circle(o1, -r1, o2, r2, &ls[2]);
  return 1;
}
// md5end

// md5begin intersect_circle_line
char intersect_circle_line(Point o, double r, Line l, Point *a1, Point *a2)
// Requires (local) vec_sub, cross, perp_line, intersect_given_x;
// (inherited) parallel_line, full_dot
{
  Point x = vec_sub(cross(perp_line(o, l), l), o);
  return intersect_given_x(o, r, x, a1, a2);
}
// md5end

/* CONVEX HULL */

// md5begin graham_cmp
int graham_cmp(const void *a, const void *b)
// Requires (local) side_of; (inherited) vec_sub, cart_cross
{
  Point p1 = *((Point *) a), p2 = *((Point *) b);
  if(p1.x == 0 && p1.y == 0)
    return -1;
  if(p2.x == 0 && p2.y == 0)
    return 1;
  ptype r = side_of(mk_point(0, 0, 1), p1, p2);
  if(r == 0)
    r = (p1.x * p1.x + p1.y * p1.y) * p2.z * p2.z
        - (p2.x * p2.x + p2.y * p2.y) * p1.z * p1.z;
  return (r < 0) ? -1 : (r > 0) ? 1 : 0;
}
// md5end

// md5begin convex_hull
void convex_hull(Point *pts, int no_points, Point *hull, int *hull_size)
// Requires (local) vec_sub, vec_add, graham_cmp; (inherited) side_of,
// cart_cross
// NOTE does not include points on the boundary (except maybe in FP)
{
  Point o;
  int ii;

  for(ii = 0, o = pts[0]; ii < no_points; ii++)
    if(pts[ii].x * o.z > pts[ii].z * o.x
        || (pts[ii].x * o.z == pts[ii].z * o.x
        && pts[ii].y * o.z > pts[ii].z * o.y))
      o = pts[ii];
  for(ii = 0; ii < no_points; ii++)
    pts[ii] = vec_sub(pts[ii], o);
  qsort(pts, no_points, sizeof(Point), graham_cmp);

  *hull_size = 0;
  for(ii = 0; ii < no_points; ii++)
  {
    while(*hull_size >= 2 &&
        side_of(hull[(*hull_size) - 2], hull[(*hull_size) - 1], pts[ii]) >= 0)
      (*hull_size)--;
    hull[*hull_size] = pts[ii];
    (*hull_size)++;
  }

  for(ii = 0; ii < *hull_size; ii++)
    hull[ii] = vec_add(hull[ii], o);
}
// md5end
